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ABSTRACT 


This  paper  Is  the  final  in  a  series  of  three  in  which  we  have  discussed 
a  finite  element  post-processing  technique.  Here  we  shall  deal  with  the 
questions  of  adaptive  mesh  selection  and  a  posteriori  error  estimation. 

Some  numerical  examples  computed  by  the  FEARS  program  will  be  used  to  Il¬ 
lustrate  the  approaches  taken. 


§1.  Introduction 


This  Is  the  final  in  a  series  of  three  papers  in  which  we  have  sought 
to  show  liow  a  suitable  post-processing  of  a  finite  element  solution  c.iii 
yield  accurate  pointwise  values  for  quantities  such  as  displacements, 
stresses,  flow  rates  and  stress  intensity  factors.  In  [1]  and  [2]  we  de¬ 
rived  a  number  of  extraction  expressions  for  such  quantities  in  the  setting 
of  some  simple  model  problems,  and  saw  how  these  expressions  could  serve  as 
the  bases  of  effective  post-processing  techniques.  We  also  carried  out 
an  error  analysis  for  such  post-processing  computations.  This  analysis: 
showed  that  the  accuracy  of  the  post-processed  value  could  be  related  to 
how  well  the  space  of  finite  element  functions  is  able  to  approximate  both 
the  solution  of  the  basic  problem  and  the  solution  of  a  related  auxiliary 
problem.  This  auxiliary  problem  is  of  the  same  form  as  the  basic  problem, 
though  with  different  loading  data.  (See  §2.5  and  §3.4  of  [1];  and  §4  of 
[2].) 

Hitherto,  except  for  a  few  qualitative  remarks,  we  have  said  litt: 
concerning  the  issues  of 

(1)  choosing  a  finite  element  subspace  for  calculating  the  approxi¬ 
mate  solution  which  is  to  be  subsequently  post-processed ; 

(ii)  estimating,  a  posteriori,  the  error  in  a  computed  post-processed 
value. 

The  significance  of  both  (1)  and  (ii)  is  quite  clear.  In  practice, 
the  goal  of  any  post-processing  computation  is  to  obtain  a  post-processed 
value  of  a  specified  accuracy  at  a  minimal  total  computational  cost.  An 
estimate  as  in  (11)  provides  a  means  of  determining  when  the  specified 
accuracy  has  been  attained,  whereas  the  choice  in  (1)  largely  determines 
the  efficiency  of  the  overall  numerical  procedure. 
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In  this  paper  we  propose  a  post-processing  algorithm.  It  is  based 
on  the  extraction  techniques  of  [1]  and  [2],  and  includes  features  that 
enable  (i)  and  (11)  to  be  handled  quite  effectively.  Our  discussion  will 
be  in  the  context  of  the  "membrane”  model  problems  already  introduced  ir  ^5 
of  [1]  and  §6  of  [2].  In  these  the  order  of  the  elements  employed  if 
fixed  (square  bilinear  elements  are  used),  and  (!)  becomes  a  matter  of 
choosing  a  finite  element  mesh  (i.e.,  nodal  points).  The  manner  in  which 
we  realize  (i)  and  (ii)  will  make  use  of  some  of  the  features  available  in 
the  FEARS  program.  In  §2  of  this  paper  we  briefly  describe  the  relevant 
features  of  FEARS.  In §3.1  we  review  the  error  analysis  of  our  earlier  papers 
[1]  and  [2]  and  show  how  this  analysis  suggests  some  a  posteriori  error 
estimates  that  can  be  computed  within  the  FEARS  framework.  These  a 
posteriori  estimates  are  the  basis  of  the  proposed  algorithm, which  we  des¬ 
cribe  in  §3.2.  Finally,  in  §4  we  return  to  some  of  the  nuroerial  examples 
of  §5  of  (1]  and  §6  of  [2],  this  time  concentrating  on  some  new  aspects 
which  are  related  to  the  Issues  (i)  and  (ii) . 


§2.  The  FEARS  program 

FEARS  is  a  research  oriented,  adaptive  finite  element  code  develojied 
at  the  University  of  Maryland.  A  detailed  description  of  the  operation 
of  the  program  can  be  found  in  [3].  For  the  purposes  of  this  paper,  the 
following  few  remarks  will  suffice.  As  already  explained  in  [2],  FEARS 
assumes  that  the  region  under  consideration  has  firstly  been  partitioned 
into  a  number  of  subregions,  each  of  which  is  a  curvilinear  quadrilateral 
Within  the  program,  each  of  these  subregions  is  transformed  by  a  change 
of  coordinates  into  a  unit  square.  The  actual  finite  element  modelling 
is  carried  out  on  these  transformed  squares.  Square  bilinear  elements 
are  used.  FEARS  has  an  adaptive  character:  starting  from  an  initial 
coarse  mesh  (usually  uniform  on  each  of  the  transformed  squares) ,  the 
program  automatically  selects,  in  a  recursive  fashion,  a  sequence  of 
"optimal"  mesh  refinements. 

The  mesh  refinement  procedure  is  based  upon  a  set  of  non-negativi 
error  indicators.  An  indicator,  say,  is  associated  with  each 

element  A.  It  is  calculated  only  using  information  about  the  finite 
element  solution  on  A  itself  and  on  the  immediately  adjacent  elements. 
These  error  indicators,  when  summed  over  all  elements  to  obtain 


(2.1) 


c  = 


L 


all  elements 

A 


say,  yield  an  estimate  for  some  user  specified  measure  of  the  error  tn 
the  finite  element  solution.  Typically,  this  measure  is  closely  related 
to  the  energy  of  the  error.  (Later  we  shall  say  a  little  more.  See  3 
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and  §4.)  Each  step  in  the  refinement  process  is  directed  towards  minimizing 
the  sum  e  of  (2.1)  in  some  "optimal"  fashion.  To  this  end,  all  elements 
A  of  an  existing  mesh  whose  error  indicators  exceed  a  threshold  value 

are  subdivided.  This  threshold  is  determined  from  some  information  on  Che 
past  history  of  the  refinement  process.  Of  course,  the  character  of  tl>e 
meshes  constructed  by  FEARS  will  depend  upon  which  error  measure  has  been 
specified.  This  reflects  the  fact  that  the  quality  of  a  mesh  is  not  an 
absolute  property,  but  must  be  viewed  relative  to  the  ultimate  goal  of  the 
finite  element  calculation.  We  shall  elaborate  further  upon  the  point  in 
§  4. 

In  FEARS  the  finite  element  equations  are  solved  by  a  direct  method. 
Usually,  after  each  refinement  step  a  full  new  solution  is  calculated  and 
a  new  set  of  error  Indicators  is  found.  However,  calculating  a  new  solu¬ 
tion  each  time  is  quite  expensive,  and  as  it  is  just  an  intermediate  step, 
only  being  used  to  compute  the  new  error  indicators  for  the  next  refine¬ 
ment  step,  one  would  like  to  avoid  it,  if  possible.  FEARS  has  a  number 
of  "economy"  modes  which  do  this  to  varying  degrees.  In  these  a  full 
solution  is  computed  only  after  a  specified  Increase  in  the  total  number 
of  elements  has  occurred  since  the  last  full  solution.  After  any  refine¬ 
ment  step  between  two  such  full  solution  phases,  the  new  error  indicators 
are  only  approximated  on  the  basis  of  the  past  history  of  the  local 
refinement  process.  This  "economy"  mode  permits  a  multi-level  refinement 
to  take  place  between  two  full  solutions.  This  possibility  is  partic¬ 
ularly  desirable  for  efficient  operation  for  problems  with  severe  singu- 

/ 

laritles.  Provided  the  number  of  such  "short  passes"  between  full  solu¬ 
tion  steps  is  not  too  great,  the  resulting  refinement  pattern  does  not 


differ  too  much  from  that  obtained  by  the  all  full  solution  method,  but, 
of  course,  with  a  considerable  saving  in  computalonal  cost. 

Estimates  such  as  e  of  (2.1),  as  well  as  being  the  basis  of  the 
automatic  mesh  refinement  feature  of  FEARS,  also  provide  a  means  of  a 
posteriori  error  estimation.  Of  course,  we  cannot  expect  that  the  esti¬ 
mate  e  should  yield  the  exact  value  of  the  specified  error  measure. 
However,  under  suitable  assumptions, estimates  may  be  calculated  that  are 
asymptotically  exact .  That  is,  if  e  denotes  the  exact  value  of  the  de¬ 
sired  error  measure,  then 

e  =  e(l+o(l))  as  e  ->■  0. 

To  illustrate  some  of  the  points  made  above,  let  us  return  to 
Example  A  of  §6  of  [2].  In  that  example  we  considered  the  boundary  value 
problem 


V^o) 

=  0 

in 

U) 

=  0 

on 

9u) 

3n 

=  0 

on 

^2 

3u) 

9n 

=  ^^2 

on 

^3’ 

where  Q  is  the  unit  circle  slit  along  the  positive  x^  axis  and 

and  are  as  shown  in  Fig.  1.  In  [2]  we  reported  on  a  sequence 

f'f  five  adaptively  refined  meshes  that  FEARS  constructed  for  this  problem 
(see  §6.2  of  [2]).  We  can  now  be  more  specific.  The  mesh  refinement 
process  for  this  example  relied  upon  error  indicators  n,  for  which  the 
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Figure  1 

The  region  for  the  model  problem  (2.2) 


sum  E  of  (2.1)  is  an  asymptotically  exact  estimate  for  the  strain  energy 

2 

measure  e  =  |v(a)-w)|  dA  of  the  error  in  the  finite  element  solution 

46 

cj.  As  we  saw  in  [2],  the  meshes  so  constructed  yield  a  rate  of  conver- 

1/2 

gence  for  the  strain  energy  norm  e  of  the  error  that  is  close  to  the 

-1/2 

theoretically  optimal  rate  of  0(N  ) ,  where  N  is  the  number  of 

degrees-of-freedom  of  the  finite  element  model.  In  contrast,  uniform 

—1/8 

meshes  would  only  give  an  0(N  )  rate.  In  Table  1  we  have  listed 

1/2  1/2 

e  and  r  for  each  of  our  meshes.  Notice  the 

e 

converging  to  1,  as  it  should,  since  e  is  an  asymptotically  exact 
estimate  for  e.  For  this  problem  FEARS  was  executed  in  an  "economy" 
mode,  so  permitting  multi-level  refinements  between  full  solutions  (the 
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§3.1.  A  Posteriori  Error  Estimates 

For  the  model  problems  discussed  in  §5  of  [1]  and  in  [2]  we  saw  that 
the  difference  between  the  exact  value  <{>  of  some  quantity  (e.g.  dis¬ 
placement,  stress,  stress  intensity  factor)  associated  with  the  exact 
solution  (jj  of  a  problem  and  an  approximate  value  i  obtained  by  suit¬ 
ably  post-processing  the  finite  element  solution  H  could  be  expressed 
as 


(3.1) 


(j,  -  $  = 


V(a)-a))  •V(iJj-iij)dA, 


where  'p  and  are  the  exact  and  finite  element  solutions  respectively 
of  an  auxiliary  problem.  This  auxiliary  problem  is  of  the  same  form  as 
the  basic  problem  for  to,  though  with  different  loading  data.  The  load¬ 
ing  data  is  determined  by  the  particular  extraction  functions  used  in  the 
post-processing  calculations. 

For  any  function  u,  let  E(u)  denote  the  "membrane"  strain  energy 
of  u. 


E(u) 


2 

|Vu|  dA. 

■  r. 


After  a  little  algebraic  manipulation  (3.1)  may  be  rewritten  as 


(3.2) 


5  -  ?  -  [E(((o-hi))  ~  (u>-hp))  -  E((ui-ip)  -  (lo-ip))]- 


Furthermore,  the  following  inequalities  also  follow  from  (3.1), 


(3.3a)  I  —  I’ I  •'  E((aj— lo)  ^  E(t()— iji)  ^  , 

(3.3b)  Id  -  $1  <  ^  (E(a,-(:,)  +  a^E(D-D)) 


for  any  real  number  a  >  0.  Notice  that  whereas  (3.2)  Is  an  exact  ex¬ 
pression  for  the  error  In  general  (3.3a)  and  (3.3b)  only  yield 

upper  bounds  for  14>-$1.  They  fail  to  take  account  of  any  cancellation 
in  the  integral  (3. 1) .  The  Inequality  (3.3a)  becomes  an  equality  only  if 
(o)  -'jj)  and  (ip  -  ill)  are  multiples  of  one  another.  (The  inequality  (3.3b) 
is  a  equality  only  in  the  more  particular  case  (w-uO  =  ^cx(i|)-iji) .) 

One  way  to  think  of  this  cancellation  phenomena  is  in  terms  of  the 
"angle"  between  the  errors  (w-u)  and  (■+'-i|)) .  Let  y  be  the  angle 
lying  between  0®  and  90°  for  which 


(3.4) 


I V(a)-w)  •V(ifj-iii)dA 
E(uj-w)^^^E(ip-iii)^^^ 


Then,  y  =  0° (cos  y=l)  corresponds  to  the  case  for  which  (3.3a)  is  an 

equality,  while  y  =  90°  (cos  indicates  complete  cancellation  in 

the  integral  (3.1).  Continuing  the  geometrical  analogy,  we  could  say 

that  in  the  case  y  =  0°,  the  errors  (w-w)  and  (ijj-ili)  are  "parallel," 

and  in  the  case  y  =  90°  they  are  "perpendicular." 

Means  of  estimating  the  quantities  appearing  on  the  right  hand  sides 

of  (3.2)  and  (3.3)  are  available  in  FEARS.  If  u  is  the  finite  element 

solution  of  a  problem  whose  exact  solution  is  u,  let  e^(u)  =  Y  n?(u) 

A  ‘ 

be  an  asymptotically  exact  estimate  for  E(u-u)  which  is  constructed 
from  elementary  error  indicators  n^(u).  That  is. 


E(u-u) 


c°(u) (1+0(1)) 


as  E(u-u)  ->  0 
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(3.5)  $  -  $  =  ((.)+;)  -  E^Oii  v)]  +  o  (1)  +  c^(ui-i/))], 


and  (3.3)  leads  to 

(3.6a)  1$  -  ?j  <  e^(<:.)^^^  e°(ij;)^^^  (l+o(l)) 

(3.6b)  \i  -  i\  <  ^(c^(w)  +  a^e°(^))  (1+0(1)) 


where  the  o(l)  term  is  valid  as  E(ai-(j)  +  E(ijj-i|i)  ->  0. 
The  equation  (3.5)  suggests  that  the  quantity 


(3.7) 


^  [G^(a)+ii^)  -  e^(u)-i|;)] 


should  provide  a  good  estimate  for  the  error  $  -  $  provided  the  o(l) 
term  is  negligible.  However,  if  the  o (1)  term  in  (3.5)  is  comparable 
•T[e^(.j+iij)-e^(u-iti)  ] 

to  M  =  — pj - pj -  then  e,  is  no  longer  reliable.  In  cases 

G  (a)+iij)+e  (oj-iji) 

where  the  angle  y  defined  in  (3.4)  is  close  to  90®,  then  u  may  be  quite 
small  and  Gj^  could  well  perform  poorly.  Turning  now  to  (3.6),  we  see 
that,  at  least  asymptotically,  the  quantities 


(3.8a) 


(3.8b)  ^3  ~  ^  (e^(d)) +a^G*^(iJi)) 

give  upper  bounds  for  |  $  -  3>|.  From  what  we  have  said  before  it  can  be 
seen  that  the  extent  of  the  asymptotic  overestimation  of  1$  -  1)j  by 
is  closely  related  to  the  angle  y  ^'or  y  near  0°,  ^2  i® 


The  Numerical  Algorithm 


Let  us  now  briefly  list  some  features  of  the  algorithm  that  we  men¬ 
tioned  in  the  Introduction: 

(a)  Instead  of  only  solving  for  w,  solutions  for  both  di  and  i]} 
are  calculated.  This  is  not  as  labourious  as  it  may  at  first  seem,  since 
the  auxiliary  problem  for  ijj  and  the  basic  problem  for  w  differ  only 
in  their  loading  data.  Thus,  they  may  be  thought  of  as  the  solutions  of 
a  multiple  load  problem  . 

(b)  The  pair  of  solutions  di  and  ij;  are  computed  for  a  sequence 
of  adaptively  refined  meshes  using  the  FEARS  program.  The  mesh  refine¬ 
ment  steps  are  based  on  the  particular  choice  of  error  indicators 

(3.9)  (.^)) 


where  a  >  0  is  some  user  chosen  constant.  The  choice  (3.9)  is  obviously 
suggested  by  (3.6b).  The  logic  being  that  with  this  choice  the  adaptive 
mesh  refinement  process  is  directed  by  an  "optimal"  minimization  of  the 
estimate  e^.  Recall  that  asymptotically  is,  in  general,  only  an 

upper  bound  for  j$  -  $j,  the  degree  of  overestimation  being  determined 
by,  amongst  other  factors,  the  angle  y  and  the  value  of  a.  It  would 
seem  preferable  to  employ  error  indicators  whose  sum  provided  a  sharper 
estimate  than  (3.6b).  Unfortunately,  the  two  other  estimates  and 

£2  that  we  have  available  cannot  be  expressed  in  the  form  (2.1).  In 
the  case  of  choosing  =  ;^(n^(a)'Hi')  “  n^((o-i|;))  would  not  always 

ensure  that  ^  0.  Tha  estimate  €2  cannot  be  conveniently  written  as 
an  element  by  element  sum,  though  by  a  proper  choice  of  a  the  values 


of  ^2  made  close.  We  shall  say  a  little  more  about  this 

later. 

(c)  Although  and  £2  cannot  be  used  in  the  role  described  in 

(b) ,  they  can  still  nonetheless  be  computed  as  global  quantities  and  be 
employed  as  a  posteriori  error  estimates.  As  such,  the  provide  a  means 
of  stopping  the  mesh  refinement  process  once  sufficient  accuracy  has  been 
attained.  Of  course,  would  provide  a  superior  estimate  to  C2 

long  as  the  above-mentioned  loss  of  reliability  associated  with  angles  "i 
near  90“  does  not  occur. 
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One  of  our  goals  in  that  example  was  to  find  approximate  values  for  the 
"stress"  '^2  =  at  P(1,0)  (see  Fig.  2).  In  [1]  we  presented  three 

different  extraction  expressions  for  this  stress.  These  extraction  ex¬ 
pressions  differed  in  the  way  that  they  handled  the  boundary  conditions 
for  the  generating  function.  Here  we  shall  only  further  discuss  the 
case  (c)  of  §5.3  of  [1];  that  is,  we  consider  post-processing  calcula¬ 
tions  based  on  the  extraction  expression 


(4.2) 


2 

V  ,'-.a)dA 


+  f  . 


dA 


where 


0  = 


(Xj^-1) 


L(Xj^-1)^+X2 


.  ,  Hill .  viV 

UXj-l)^+l  4+X2  ^  (. 


The  auxiliary  problem  associated  with  (4.2)  is 

2  2 

V  fp  =  -V  p  in  n 

(4.3) 

P  =0  on  sn, 

which,  as  it  should  be,  is  of  the  same  form  as  (4.1)  but  with  different 
right  hand  sides. 

In  [1]  we  saw  that  (4.2)  led  to  highly  accurate  approximations  for 
even  in  the  case  of  quite  coarse  meshes.  The  numerical  results  we 
reported  there  were  for  a  sequence  of  uniform  meshes.  At  the  time  we 
did  not  comment  on  this  particular  selection  of  meshes.  Qualitatively 
however,  such  a  choice  would  not  seem  unreasonable  from  the  viewpoint  of 


the  theory  we  developed  in  [1]  and  [2].  The  solutions  of  both  the  basic 
problem  (A.l)  and  the  auxiliary  problem  (4.3)  are  relatively  smt>o(li. 

Thus,  uniform  meshes  would  seem  appropriate.  In  fact,  the  adaptive  mesh 
refinement  algorithm  that  we  outlined  in  §3.2  exactly  reproduces  this 
sequence  of  meshes,  independently  of  n. 

Using  the  estimates  given  in  §3.1  we  are  able  to  estimate  the  error 
in  the  post-processed  value  4'^.  For  the  above  meshes.  Table  2  lists  a 
selection  of  a  posteriori  estimates  for  1^2  -  $2 i  based  on  (3.7)  and 
(3.8).  Notice  that  in  the  case  of  the  ratio  of  the  estimated  error 

to  the  true  error  in  appears  to  converge  to  1  as  the  mesh  size  is 

decreased.  This  is  consistent  with  (3.5).  On  the  other  hand,  for  i.., 
and  [. 2  the  corresponding  ratios  each  seem  to  stabilize  around  value.s 
greater  than  1  as  the  mesh  is  refined.  Again,  this  is  expected  on  the 
basis  of  (3.6)  which  shows  that  ^2  and  will,  in  general,  asymptote 

to  upper  bounds  for  the  error.  Let  us  also  remark  that  since  £2  gives 
only  a  slight  overestlmatlon  of  |$2  ~  ’^’2!  behaves  quite  well, 

we  should  expect  that  the  angle  y  between  (u-w)  and  (il^-ij))  is  not 
too  close  to  90°.  In  Table  2  we  have  listed  an  estimate  for  this  angle 
which  confirms  this  expectation.  Notice  also  the  Interesting  fact  that 
y  is  almost  Independent  of  the  mesh. 


No.  of  elements  in  quarter 
segment  (uniform  mesh) 


relative  error  in  standard 
finite  element  stress  value 

29% 

1 

1*2  - 

error  in  post-processed 
stress  value 

10.495(-3) 

(|<1)2  - 

(1.5%) 

Y 

angle  between  (cj-oa)  and 
(see  (3.4)) 

37.1“ 

e  , 

estimated  error 

-  $2!) 

10.351(-3) 

(.986) 

"2 

12.976(-3) 

(1.236) 

1 

^3 

24.412(-3) 

(2.326) 

2.553(-3) 

(.37%) 


2.599(-3) 

(1.018) 

3.255(-3) 

(1.275) 

6.444(-3) 

(2.524) 


0.637(-3) 


0.641(-3) 

(1.006) 

0.820(-3) 

(1.287) 

1.689(-3) 

(2.651) 


TABLE  2.  A  posteriori  estimates  of  the  error  in  I  ^  for  the  example  lis 
cussed  in  §4.1.  =  ■^|c^((ii+^)  -  c^((ji-'^)  j  ,  f.  ^  (;>) 

e,  =  4(e°(a')) +e°(ii)). 


1|  0.~  .7..  Of-  ,  0.^.  1/2  0,  . 

y\f  (ui+\Ij)  -  c  (lii-v)  I  ,  =  f  (I'O  €  (  ,') 


§4.2.  The  Slit  Membrane  Example  of  §6.2  of  [2] 

In  this  section  we  shall  return  to  Example  A  of  §6  of  [2).  We  have 
already  used  this  example  in  g 2  of  this  paper  to  illustrate  some  of  the 
features  of  the  FEARS  program.  The  governing  equations  and  boundary  con¬ 
ditions  are  given  in  (2.2)  and  depicted  in  Fig.  1.  As  in  [2],  let  us  be 
Interested  in  using  post-processing  techniques  to  find  approximations  to 
the  leading  stress  intensity  factor  kj .  Recall  that  was  defined  to 

be  the  coefficient  of  the  leading  term  in  the  asymptotic  expansion 

(4.4)  M  =  k  r^^^  sin  y  +  (k  =  1.35812) 

14  1 


for  u>  in  the  vicinity  of  the  slit  tip  (0,0).  Here  we  shall  only  treat 
post-processing  calculations  based  upon  the  extraction  expression 


(4.5) 


P  -  '7'*'  ’tidOds 


where 


2  -1/4 

TI 


The  expression  (4.5)  is  an  instance  of  what  we  have  been  calling  a 
generalized  influence  function  method.  The  auxiliary  function  Introduced 
in  the  error  analysis  of  the  extraction  expression  (4.5)  is 
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By  chance.  In  this  case,  we  can  explicitly  solve  (4.6)  to  obtain 


(4.7) 


,  21/4.0 

'!)  =  —  r  sin  T  , 

T  4 


which  is  just  a  multiple  of  the  leading  term  in  the  expansion  (4.4)  of 
u).  (Of  course,  we  cannot  expect  to  be  able  to  do  this  so  simply  in 
general . ) 

As  explained  in  §2,  the  meshes  we  considered  there  and  in  [2]  for 
this  problem  were  constructed  by  FEARS  using  the  error  indicator  = 
n^(iij)  •  (That  is,  on  the  basis  of  minimizing  the  estimate  of  the  strain 

energy  E(to-di)  of  the  error.)  Strictly  speaking,  this  indicator  does  not 
fit  into  the  framework  we  outlined  in  (b)  of  §3.2;  though,  it  can  be 
thought  of  as  a  limiting  case  as  a  •>  0.  However,  executing  the  algorithm 
of  §3.2  with  a  number  of  choices  of  a  leads  to  sequences  of  meshes  which, 
though  different  from  those  above,  show  much  the  same  refinement  charac¬ 
teristics.  For  this  reason,  and  since  we  want  our  numerical  results  here 
to  complement  those  of  [2],  we  shall  work  with  the  same  sequence  of  meshes 
constructed  for  this  problem  in  [2]. 

To  try  to  see  why  the  character  of  the  meshes  constructed  above  should 
be  Independent  of  a  let  us  rewrite  (4.4)  as 


3/4 

where  ojq  =  0(r  ).  Let  (Oq  be  the  finite  element  solution  that  would 
be  obtained  were  the  loading  in  (2.2)  such  that  Uq  was  the  exact  solu¬ 
tion.  By  the  linearity  of  the  model  problem. 


The  function  is  relatively  smooth  in  comparison  with  ip,  however, 

from  our  point  of  view  the  relative  magnitude  of  the  factor  ^  is 
equally  important.  If  it  is  sufficiently  large,  then  ('p-i)  will 

make  the  major  contribution  to  the  error  to  -  w;  while  if  it  is  small 


enough,  then  for  meshes  that  are  not  too  fine,  the  to^  -  uJq  contribution 
will  dominate.  As  long  as  the  first  case  applies,  then  an  "ideal"  finite 
element  mesh  would  exhibit  a  severe  refinement  about  the  slit  tip  at 
(0,0),  while  if  the  latter  case  applies  a  more  uniform  mesh  is  called 
for,  at  least  to  begin  with.  Which  case  actually  occurs  in  this  problem 
is  clearly  demonstrated  by  the  results  shown  in  Table  3 — for  the  level  of 
refinement  encountered  there  the  -  to^  contribution  is  relatively 
small.  It  is  not  surprising  then,  that  regardless  of  whether  the  mesh 

refinement  process  is  directed  towards  minimizing  ECw-w),  or 

~  2  ' 

some  combination  ECui-w)  +  a  E(\p~<p) ,  the  resulting  meshes  will  not  be 
significantly  different.  However,  we  should  point  out,  and  the  trend  in 
Table  3  shows  this,  that  at  some  point  the  magnitude  of  the  two  components 
in  (4.9)  will  become  comparable.  At  that  point,  the  different  choices  of 
a  will  lead  to  meshes  of  significantly  different  natures. 
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Mesh  Label 

|k^  -1  1 

1 

E(a>— to)  ^ 

I 

00 

.21 

II 

1.02 

.22  i 

} 

III 

.99 

.25  1 

IV 

.96 

.34 

V 

.94 

.38 

1 

TABLE  3.  Magnitude  of  the  ~  (4'-'}')  and  components 

of  the  error  (w-tZ) . 

In  Table  4  we  list  a  number  of  different  a  posteriori  estimates  ot 
the  error  in  for  the  above  sequence  of  meshes,  along  with  the  exact 

values  of  the  error  and  the  angle  y  which  are  able  to  be  found  analyt¬ 
ically  for  this  particular  problem.  The  apparent  asymptotic  exactness 
of  is  consistent  with  (3.5)  and  the  fact  that  the  angle  y  is 

rather  small.  The  estimates  and  appear  to  lead  to  upper  bounds 

for  |kj^-kj^l.  Since  the  angle  y  is  near  0®  we  would  expect  that  c  ^ 
would,  in  the  limit,  give  only  a  slight  overestimation  of  the  error. 

Again  the  numerical  results  in  Table  4  are  consistent  with  this. 


Mesh  Label 
(degrees-of- freedom) 


1  1  1 
(118) 


IV 

(171) 


(391) 


error  In  post-processed  I 
stress  intensity  factor  I  .193860 


(|  /|  k^l  ) 


angle  between  (oj-  '.  ) 
and  (ip-tp)  (see  (3.4)) 


— 1 -  ^ 

i  i 

;  1 

:  1 

.029140  ,  .010690  ! 

(2.1%)  (.79%) 

;  i 

20.1°  :  22,4°  : 

1 

i  ! 

1 
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54.3.  A  Modification  of  the  Slit  Membrane  Example  of  ^6.2  of  [2] 

Let  us  now  alter  the  basic  problem  (2.2)  that  we  have  been  consider¬ 
ing  in  §4.2  by  modifying  the  loading  so  as  to  reduce  the  leading  stress 
intensity  factor  by  an  order  of  magnitude  or  so.  Specifically,  in  place 
of  (2.2)  consider 

2  4, 

V  (0  =  0  in  0 

w*  =  0  on  I’j 

(4.10) 


^  +'  (1.3  sin|)  on  13 

which  has  the  exact  solution  cn*  =  w  +  1.3  r^^^  sin  where  w  is  the 
solution  of  (2.2).  The  leading  stress  intensity  factor  for  co*  is 

k*  =  +  1.3  =  -.058122. 

The  algorithm  §3.2  was  executed  for  this  problem  with  the  following 
three  different  choices  of  error  indicator  of  (3.9)  governing  the 

mesh  construction  process: 

(A)  =  n^(io*)  (This  can  be  thought  of  as  a  limiting  form  of  (3.9)  as 

a  -*  0 ,  the  normalizing  factor  ~  of  (3.9)  having  no  effect  on  the  re- 

2o 

flnement  process.) 

(B)  (n°(.>)  +  i  n°('^)) 

(C)  n^  =  n^(ii').  (This  can  be  thought  as  a  limiting  form  of  (3.9)  as  a.  ->■  .) 
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Some  numerical  results  for  a  selection  of  meshes  obtained  using  each 
of  these  indicators  is  reported  in  Table  5.  The  first  thing  to  notice  from 
Table  5  is  the  very  significant  differences  between  the  meshes  created 
using  each  of  the  above  three  Indicators.  As  the  "distribution  of  elements" 
columns  show,  strategy  (C)  leads  to  meshes  that  concentrate  elements  in  the 
inner  subregions  near  the  slit  tip,  while  (B)  and  (A)  produce  progres.s  i vely 
less  severe  refinement.  This  effect  can  be  readily  explained:  In  place  of 
(4.9)  we  now  have, 

W*  -  0)*  =  ^  +  (cJq-Wq)  , 

and  since  |k*|  is  now  significantly  smaller  than  |kj|,  we  can  no  longer 
regard  the  (ojq-Uq)  contribution  to  the  error  as  negligible,  as  we  could 
in  §4.2.  Indeed,  for  the  initial  uniform  mesh  I*,  from  which  each  of  our 
strategies  (A),  (B)  and  (C)  starts, 

— - — ^ - =  .20  and  - -  ■  ^  .  =  .98 

which  is  almost  an  exact  reversal  of  the  situation  for  mesh  I  in  §4.2 
(see  Table  3).  Therefore,  it  is  to  be  expected  that  strategy  (A),  which 
seeks  to  minimize  E(w*-a)*);  strategy  (B) ,  which  seeks  to  minimize 
(E(aj*-a)*)  + -^  E(i()-iij));  and  strategy  (C)  ,  which  seeks  to  minimize 
E(i|;-ij;*)  ,  will  now  lead  to  meshes  that  are  of  a  quite  different  nature 
from  the  very  start.  The  meshes  produced  by  (C)  should  exhibit  the 
severest  refinement,  reflecting  the  fact  that  such  a  refinement  is 
necessary  to  approximate  the  singular  function  ip  well.  On  the  other 
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hand,  (A)  should  produce  the  most  uniform-llke  mesh,  since  early  on,  tlu 
smooth  component  -  Wq  dominates  the  error  ,  and  uniform-llke 

meshes  are  sufficient  to  approximate  a>Q  well.  The  actual  meshes  pro¬ 
duced  are  in  accord  with  these  expectations. 

Notice  also  that  the  meshes  produced  by  the  algorithm  seem  to  have 

achieved  their  respective  goals  in  some  "optimal"  fashion,  at  least  in 

comparison  to  one  another.  Among  the  meshes  III*,  lU?  and  III*,  each  of 

which  has  a  comparable  number  of  degrees-of-freedom,  the  mesh  IHJ  has 

the  lowest  value  for  E(t»)*-ai*)  (see  column  headed  A  ),  the  mesh  IIIJ 

A  o 

/5  ~  1 

has  the  lowest  value  for  -j-  (E((ji)*-a)*)  +—  E(ii)-iJj*))  (see  column  headed 
Ag) ,  while  the  mesh  IH^  has  the  smallest  energy  error  in  (see 

column  headed  A_) . 

Let  us  now  turn  to  the  accuracy  in  the  post-processed  value  ic*  for 
the  various  meshes.  As  we  have  said  previously,  this  accuracy  is  dependent 
not  only  on  the  accuracy  of  m*  but  also  on  the  accuracy  of  Eeferring 

to  Table  5,  we  see  that  the  meshes  II*  and  though  yielding  good 

accuracies  for  w*  in  the  energy,  are  the  worst  members  from  among  the 
sets  II*  and  III*  respectively  as  far  as  the  accuracy  of  ip*  in  energy  is 
concerned.  Not  surprisingly  then,  the  accuracy  of  ic*  for  the  meshes 
II*  and  IIIJ  is  inferior  to  that  for  the  corresponding  meshes  constructed 
using  strategies  (B)  and  (C) .  The  question  obviously  arises  as  to  what 
is  the  choice  of  a  for  use  in  (3.9)  that  gives  the  best  accuracy  in  k* 
for  a  given  number  of  degrees-of-freedom.  This  is  a  difficult  matter  to 
analyse  in  detail,  though  the  results  of  Table  5  seem  to  indicate  that, 
at  least  fr  this  problem,  the  results  are  rather  insensitive  to  non-rero 
choices  of  a.  As  a  rough  guide  though,  and  this  is  how  we  chose  our 


a=  —  ,  a  could  be  selected  so  that  for  the  Initial  uniform  mesh, 

/5 

and  a  are  about  equal.  The  logic  I'Cliind  this  elioiec 

being  that  this  a  approximately  equalizes  *^he  asymptotical  1  v 

sharper  estimate  least  for  the  initial  uniform  mesh) .  Although 

we  had  no  need  to  do  so  in  our  calculations,  one  could  possibly  change  a 
during  the  course  of  the  calculations  to  ensure  that  t'2  ^2  always 

remain  close  to  one  another. 

Finally,  let  us  say  a  little  about  the  a  posteriori  estimates  for 
|k*-kj|  based  upon  f'2  ^3"  Some  numerical  results  for  the  mesh 

sequence  I*,  II*  and  III*  are  reported  in  Table  6.  We  also  list  the 

0  D 

value  of  the  angle  y  between  the  errors  ((ii*-^*)  and  (ip-ij,).  For 
this  problem  we  are  able  to  calculate  y  exactly.  Notice  that  these 
angles  are  rather  close  to  the  critical  value  y  =  90“ .  As  we  explained 
in  §3.1,  for  such  values  of  y  we  expect  the  estimate  to  be  un¬ 

reliable  due  to  the  significance  of  the  second  term  of  (3.5).  In  addi¬ 
tion,  the  estimate  should  be  a  considerable  overestimation  of 

|k*-k*|.  Our  numerical  results  seem  to  confirm  both  these  points. 

In  general  of  course,  we  cannot  calculate  y,  and  so  in  a  case  such 

as  the  above  we  could,  just  on  the  basis  of  be  deceived  into  believ¬ 

ing  that  there  is  greater  accuracy  in  k*  than  is  actually  the  case. 
However,  the  fact  that  and  differ  by  an  order  of  magnitude  or 

so  should  act  as  a  warning  that  we  are  more  than  likely  in  a  near  critical 

^1 

case.  (In  fact,  the  ratio  —  is  an  estimate  for  cos  y.)  In  such  a  sl- 
tuatlon,  should  not  be  trusted,  while  could  be  used  to  gauge  the 

accuracy  of  fe*,  realizing  however,  that  it  is  probably  a  considerable 
overestimation  of  the  error. 


Mesh  Label 
(degrees*of- freedom) 


error  In  post-processed 
stress  intensity  factor 

|{k*-k*|/|k*i) 


.8357(-2) 

(14.4%) 


.2239(-2) 

(3.85%) 


! 

i 

angle  between  (ai*-tL*) 
and  (iii-il,)  (see  (3.4)) 

— 

78.5° 

r 

84.0° 

j 

1 

82.8° 

! 

£ 

3 

estimated  error 

(e^/|k*-k*l ) 

.2360(-2) 

.1566(-2) 

1 

.0l03(-2) 

(.28) 

(.70) 

(.14) 

"2 

1.7513(-2) 

1.4005(-2) 

.5170(-2) 

(2.10) 

(6.26) 

(7.21) 

"3 

1.9l50(-2) 

1.4042(-2) 

.5174(-2) 

1 

(2.29) 

(6.27) 

(7.22'' 

TABLE  6.  A  posteriori  estimates  for  jkj-kj)  for  the  example  of  §4.3. 
£  =  ^1  C  ^  (ai  )  -  f  ^{uo*-ii))  1  ,  =  £*^(u'*)^^^  £^(lii). 


£3  =  ^  (cO(uA)  +  \  z^a)). 
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§ 5.  Concluding  Remarks 

Our  main  intention  in  this  paper  was  to  examine  the  behaviour  of  the 
algorithm  of  §3.2  in  some  model  post-processing  applications,  concentrating 
especially  on  the  mesh  selection  and  error  estimation  features  of  the  algorithm. 
As  far  as  accuracy  of  the  post-processed  value  was  concerned,  the  particular 
examples  dealt  with  here  proved  rather  insensitive  to  the  choice 
of  a  in  the  error  indicator  (3.9),  In  our  most  extreme  example  the  spread 
in  accuracy  for  the  same  number  of  degrees-of-f reedom,  was  only  by  a  fac¬ 
tor  of  2.  Whether  this  Insensitivity  is  a  property  to  be  expected  in 
general  is  an  open  question  at  the  moment.  However,  we  can  at  least  say 
from  a  theoretical  point  of  view  that  some  dependence  on  a  is  to  be 
expected  whenever  the  basic  and  auxiliary  problems  have  solutions  with  dif¬ 
ferent  "smoothness"  characteristics. 

Our  examples  also  show  that,  except  for  the  critical  casp  when  ^  is 
near  90°,  the  error  estimates  and  perform  well  asymptoticaJ  ly . 

The  occurrence  of  the  critical  case  can  be  detected  numerically  by  the 
^2 

fact  that  —  ^  5,  say.  In  this  critical  case,  even  though  c  appears 
^1 

unreliable,  still  seems  to  provide  a  useable  estimate,  albeit  an 


overly  pessimistic  one. 
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•  To  help  bridge  gaps  between  computational  directions  in  engineering, 
physics,  etc.  and  those  in  the  mathematical  community. 

•  To  provide  a  limited  consulting  service  in  all  areas  of  numerical 
mathematics  to  the  University  as  a  whole,  and  also  to  government 
agencies  and  industries  in  the  State  of  Maryland  and  the  Washington 
Metropolitan  area. 

•  To  assist  with  the  education  of  numerical  analysts,  especially  at  the 
postdoctoral  level,  in  conjunction  with  the  Interdisciplinary  Applied 
Mathematics  Program  and  the  programs  of  the  Mathematics  and  Computer 
Science  Departments.  This  Includes  active  collaboration  with  government 
agencies  such  as  the  National  Bureau  of  Standards. 

•  To  be  an  International  center  of  study  and  research  for  foreign 
students  In  numerical  mathematics  who  are  supported  by  foreign 
governments  or  exchange  agencies  (Fulbrlght,  etc.). 

Further  Information  may  be  obtained  from  Professor  I.  BabuSka,  Chairman, 
Laboratory  for  Numerical  Analysis,  Institute  for  Physical  Science  and 
Technology,  University  of  Maryland,  College  Park,  Maryland  20742. 


